#########################################
## Cell Med Review                     ##
## Gastric Bypass/ UOS                 ##  
## MLM: Diabetes (HOMAIR/B) and Weight ##
## Ben Krick                           ## 
## January 27th, 2022                  ##
#########################################
library(lme4)
library(lmerTest)
library(sjPlot)
library(jtools)
library(ggplot2)
library(readr)
library(ggeffects)
library(broom)
library(tidyverse)

setwd("/Volumes/dfs/Playdon Research Group/Team_Members_Students/Ben Krick/UOS: GB : Ceramide Study /3_Projects/Gastric_Bypass_Cell_Med")

# load in the data 
GB_data_long <- read_csv("Data/GB_data_long.csv")
#GB_data_long <- read_csv("GB_data_long_july_28.csv")
GB_data_long$group <- as.factor(GB_data_long$group)


## identify the ceramides you are using for analysis 
ceramides <- c("Cer(d18:0/16:0)","Cer(d18:0/18:0)", "Cer(d18:0/20:0)", "Cer(d18:0/22:0)","Cer(d18:0/24:0)",
               "Cer(d18:0/24:1)", "Cer(d18:1/16:0)", "Cer(d18:1/18:0)", "Cer(d18:1/20:0)", "Cer(d18:1/22:0)",
               "Cer(d18:1/24:0)", "Cer(d18:1/24:1)")

#log and save names 
GB_data_long[198:209] <- lapply(GB_data_long[ceramides], function (x) ifelse(x != 0, log10(x), x))

colnames(GB_data_long)[198:209] <- paste("log", colnames(GB_data_long)[198:209],sep = "_")

log_ceramides <- colnames(GB_data_long)[198:209]

# # # Mixed Effect Model: Ceramides and hba1c

#create empty final results 
results <- data.frame(matrix(ncol = 4, nrow = 0))

# run the model in a loop of ceramides (original model in paper)
for (i in log_ceramides){
  exposure = i 
  model <- lmer(log10(hba1c) ~ get(exposure) + group +  bmi + age  + education + sexint + marrstat + ethnicity + income  + bpmeds + bcholmed + triglyceride + cholesterol 
                + (time|unqid), REML=F, 
                na.action = na.exclude, 
                data=GB_data_long)
  temp_results<- cbind(summary(model)[["coefficients"]][2,1],
                       confint.merMod(model, method = "Wald")[6,1],
                       confint.merMod(model, method = "Wald")[6,2],
                       summary(model)[["coefficients"]][2,5])
  row.names(temp_results) <- i 
  results <- rbind(results, temp_results)
}


#provide column names
colnames(results) <- c("ESTIMATE", "CI_L", "CI_U", "PVALUE")
#add BH and organize 
results$BH <- p.adjust(results$PVALUE,  method = "BH") #adjusted p value 

results <- results %>% 
  rownames_to_column(var="ceramide")%>% 
  remove_rownames

write_csv(results, "Results/supp_table5a.csv")

#loop through the figures 
# dat1 <- ggpredict(model1, terms = "log10_cer18_16")
# plot(dat1)

#create empty final results 
results <- data.frame(matrix(ncol = 4, nrow = 0))

# run the model again in a loop of ceramides (just sugery group, remove control for group)
for (i in log_ceramides){
  exposure = i 
  model <- lmer(log10(hba1c) ~ get(exposure) +  bmi + age  + education + sexint + marrstat + ethnicity + income  + bpmeds + bcholmed + triglyceride + cholesterol 
                + (time|unqid), REML=F, 
                na.action = na.exclude, 
                data=GB_data_long[GB_data_long$group %in% "surgery",])
  temp_results<- cbind(summary(model)[["coefficients"]][2,1],
                       confint.merMod(model, method = "Wald")[6,1],
                       confint.merMod(model, method = "Wald")[6,2],
                       summary(model)[["coefficients"]][2,5])
  row.names(temp_results) <- i 
  results <- rbind(results, temp_results)
}


#provide column names
colnames(results) <- c("ESTIMATE", "CI_L", "CI_U", "PVALUE")
#add BH and organize 
results$BH <- p.adjust(results$PVALUE,  method = "BH") #adjusted p value 

# # # results are not significant -- maybe because of power or is nonsurgery driving the relationship
results <- results %>% 
  rownames_to_column(var="ceramide")%>% 
  remove_rownames



# # # Mixed Effect Model: Ceramides and homa_IR


# run the model in a loop of ceramides 

#create empty final results 
results <- data.frame(matrix(ncol = 4, nrow = 0))


# run the model in a loop of ceramides (original model in paper)
for (i in log_ceramides){
  exposure = i 
  model <- lmer(log10(homair) ~ get(exposure) + group +  bmi + age  + education + sexint + marrstat + ethnicity + income  + bpmeds + bcholmed + triglyceride + cholesterol 
                + (time|unqid), REML=F, 
                na.action = na.exclude, 
                data=GB_data_long)
  temp_results<- cbind(summary(model)[["coefficients"]][2,1],
                       confint.merMod(model, method = "Wald")[6,1],
                       confint.merMod(model, method = "Wald")[6,2],
                       summary(model)[["coefficients"]][2,5])
  row.names(temp_results) <- i 
  results <- rbind(results, temp_results)
}


#provide column names
colnames(results) <- c("ESTIMATE", "CI_L", "CI_U", "PVALUE")
#add BH and organize 
results$BH <- p.adjust(results$PVALUE,  method = "BH") #adjusted p value 

results <- results %>% 
  rownames_to_column(var="ceramide")%>% 
  remove_rownames

write_csv(results, "Results/supp_table5b.csv")



#create empty final results 
results <- data.frame(matrix(ncol = 4, nrow = 0))

# run the model again in a loop of ceramides (just sugery group, remove control for group)
for (i in log_ceramides){
  exposure = i 
  model <- lmer(log10(homair) ~ get(exposure) +  bmi + age  + education + sexint + marrstat + ethnicity + income  + bpmeds + bcholmed + triglyceride + cholesterol 
                + (time|unqid), REML=F, 
                na.action = na.exclude, 
                data=GB_data_long[GB_data_long$group %in% "surgery",])
  temp_results<- cbind(summary(model)[["coefficients"]][2,1],
                       confint.merMod(model, method = "Wald")[6,1],
                       confint.merMod(model, method = "Wald")[6,2],
                       summary(model)[["coefficients"]][2,5])
  row.names(temp_results) <- i 
  results <- rbind(results, temp_results)
}


#provide column names
colnames(results) <- c("ESTIMATE", "CI_L", "CI_U", "PVALUE")
#add BH and organize 
results$BH <- p.adjust(results$PVALUE,  method = "BH") #adjusted p value 

results <- results %>% 
  rownames_to_column(var="ceramide")%>% 
  remove_rownames



# # # Mixed Effect Model: Ceramides and HOMA-B

# # # will have to add this in 
homab <- read_csv("/Volumes/dfs/Playdon Research Group/Team_Members_Students/Ben Krick/UOS: GB : Ceramide Study /1_Raw Data /UOS_covariates.csv")

time1 <- homab[c("unqid","homab_1")]
colnames(time1)[2] <- "homab"
time1$time <- 0 

time2 <- homab[c("unqid","homab_2")]
colnames(time2)[2] <- "homab"
time2$time <- 2

time4 <- homab[c("unqid","homab_4")]
colnames(time4)[2] <- "homab"
time4$time <- 12

homab <- rbind(time1, time2, time4)

GB_data_long <- merge(GB_data_long, homab)


#create empty final results 
results <- data.frame(matrix(ncol = 4, nrow = 0))


# run the model in a loop of ceramides (original model in paper)
for (i in log_ceramides){
  exposure = i 
  model <- lmer(log10(homab) ~ get(exposure) + group +  bmi + age  + education + sexint + marrstat + ethnicity + income  + bpmeds + bcholmed + triglyceride + cholesterol 
                + (time|unqid), REML=F, 
                na.action = na.exclude, 
                data=GB_data_long)
  temp_results<- cbind(summary(model)[["coefficients"]][2,1],
                       confint.merMod(model, method = "Wald")[6,1],
                       confint.merMod(model, method = "Wald")[6,2],
                       summary(model)[["coefficients"]][2,5])
  row.names(temp_results) <- i 
  results <- rbind(results, temp_results)
}


#provide column names
colnames(results) <- c("ESTIMATE", "CI_L", "CI_U", "PVALUE")
#add BH and organize 
results$BH <- p.adjust(results$PVALUE,  method = "BH") #adjusted p value 

results <- results %>% 
  rownames_to_column(var="ceramide")%>% 
  remove_rownames

write_csv(results, "Results/ceramides_homab.csv")






# instead of running for bmi (change at each time points)
